% Producing Figure 10

close all
warning off

load('no_policy.mat')

n1  = run2_avg_n(1);
k1  = run2_avg_k(1);
rb1 = run2_avg_rb(1);

t_before = run_before;
t_after  = run_after;

tfp_shock = run2_avg_tfp;

eta1_rhsee_nopol = solved_eta1_rhsee;
eta2_rhsee_nopol = solved_eta2_rhsee;
eta3_rhsee_nopol = solved_eta3_rhsee;
eta1_rbar_nopol  = solved_eta1_rbar;
eta2_rbar_nopol  = solved_eta2_rbar;
eta3_rbar_nopol  = solved_eta3_rbar;

filename = 'policy_threL_150.mat';

temp = load(filename,'solved_eta1_rhsee');
eta1_rhsee_pol = temp.solved_eta1_rhsee;
temp = load(filename,'solved_eta2_rhsee');
eta2_rhsee_pol = temp.solved_eta2_rhsee;
temp = load(filename,'solved_eta3_rhsee');
eta3_rhsee_pol = temp.solved_eta3_rhsee;
temp = load(filename,'solved_eta4_rhsee');
eta4_rhsee_pol = temp.solved_eta4_rhsee;

temp = load(filename,'solved_eta1_rbar');
eta1_rbar_pol = temp.solved_eta1_rbar;
temp = load(filename,'solved_eta2_rbar');
eta2_rbar_pol = temp.solved_eta2_rbar;
temp = load(filename,'solved_eta3_rbar');
eta3_rbar_pol = temp.solved_eta3_rbar;

temp   = load(filename,'thre_l');
thre_l = temp.thre_l;

filename = 'policy_tau_5_always.mat';

temp = load(filename,'solved_eta1_rhsee');
eta1_rhsee_pol3 = temp.solved_eta1_rhsee;
temp = load(filename,'solved_eta2_rhsee');
eta2_rhsee_pol3 = temp.solved_eta2_rhsee;
temp = load(filename,'solved_eta3_rhsee');
eta3_rhsee_pol3 = temp.solved_eta3_rhsee;

temp = load(filename,'solved_eta1_rbar');
eta1_rbar_pol3 = temp.solved_eta1_rbar;
temp = load(filename,'solved_eta2_rbar');
eta2_rbar_pol3 = temp.solved_eta2_rbar;
temp = load(filename,'solved_eta3_rbar');
eta3_rbar_pol3 = temp.solved_eta3_rbar;

temp = load(filename,'tau');
tau  = temp.tau;

%% simulate

options = optimset;

av   = zeros(t_before+1+t_after,1);

kv   = zeros(t_before+1+t_after,1);
rbv  = zeros(t_before+1+t_after,1);
nv   = zeros(t_before+1+t_after,1);
xiv  = zeros(t_before+1+t_after,1);
hv   = zeros(t_before+1+t_after,1);
cv   = zeros(t_before+1+t_after,1);
yv   = zeros(t_before+1+t_after,1);
iv   = zeros(t_before+1+t_after,1);
lv   = zeros(t_before+1+t_after,1);
pibv = zeros(t_before+1+t_after,1);
runv = zeros(t_before+1+t_after,1);
rhseev = zeros(t_before+1+t_after,1);
rkhat  = zeros(t_before+1+t_after,1);
rkhats = zeros(t_before+1+t_after,1);
prv  = zeros(t_before+1+t_after,1);

kv2   = zeros(t_before+1+t_after,1);
rbv2  = zeros(t_before+1+t_after,1);
nv2   = zeros(t_before+1+t_after,1);
xiv2  = zeros(t_before+1+t_after,1);
hv2   = zeros(t_before+1+t_after,1);
cv2   = zeros(t_before+1+t_after,1);
yv2   = zeros(t_before+1+t_after,1);
iv2   = zeros(t_before+1+t_after,1);
lv2   = zeros(t_before+1+t_after,1);
pibv2 = zeros(t_before+1+t_after,1);
runv2 = zeros(t_before+1+t_after,1);
rhseev2 = zeros(t_before+1+t_after,1);
rkhat2  = zeros(t_before+1+t_after,1);
rkhats2 = zeros(t_before+1+t_after,1);
prv2  = zeros(t_before+1+t_after,1);

kv3   = zeros(t_before+1+t_after,1);
rbv3  = zeros(t_before+1+t_after,1);
nv3   = zeros(t_before+1+t_after,1);
xiv3  = zeros(t_before+1+t_after,1);
hv3   = zeros(t_before+1+t_after,1);
cv3   = zeros(t_before+1+t_after,1);
yv3   = zeros(t_before+1+t_after,1);
iv3   = zeros(t_before+1+t_after,1);
lv3   = zeros(t_before+1+t_after,1);
pibv3 = zeros(t_before+1+t_after,1);
runv3 = zeros(t_before+1+t_after,1);
rhseev3 = zeros(t_before+1+t_after,1);
rkhat3  = zeros(t_before+1+t_after,1);
rkhats3 = zeros(t_before+1+t_after,1);
prv3  = zeros(t_before+1+t_after,1);
lv_unreg  = zeros(t_before+1+t_after,1);
kv_unreg  = zeros(t_before+1+t_after,1);
erk_next_unreg  = zeros(t_before+1+t_after,1);

%% Simulation

kv(1)  = k1;
rbv(1) = rb1;
nv(1)  = n1;
lv(1)  = kv(1)/nv(1);

kv2(1)  = kv(1);
rbv2(1) = rbv(1);
nv2(1)  = nv(1);
lv2(1)  = kv(1)/nv(1);

kv3(1)  = kv(1);
rbv3(1) = rbv(1);
nv3(1)  = nv(1);
lv3(1)  = kv(1)/nv(1);

    for t=2:1+t_before+1+t_after

        av(t) = tfp_shock(t-1);
        
        % no policy
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x*eta1_rhsee_nopol);
                rbv(t)    = exp(x*eta1_rbar_nopol);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x*eta3_rhsee_nopol);
                rbv(t)    = exp(x*eta3_rbar_nopol);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x*eta2_rhsee_nopol);
            rbv(t)    = exp(x*eta2_rbar_nopol);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk_next = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk_next,lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next - erk_next) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;
        
    % leverage cap policy %%%%%%%%%%%%%%%%%%
    
        hv2(t)  = ((1-alpha)/psi*exp(av(t))*kv2(t-1)^(alpha))^(1/(alpha+1/nu));
        yv2(t)  = exp(av(t)) * kv2(t-1)^(alpha)*hv2(t)^(1-alpha);

        rkhats2(t) = log ( rbv2(t-1) * (1-1/lv2(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat2(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t-1));
        xiv2(t)  = exp(rkhat2(t))/exp(rkhats2(t));

        K  = [kv2(t-1) rbv2(t-1) nv2(t-1) xiv2(t)];
        x0 = (log(K)).^expmx;
        x  = prod(x0,2)';
        if xiv2(t) >= 1
            pibv2(t)  = ( exp(rkhat2(t)) + 1-delta )*lv2(t-1)*nv2(t-1) - rbv2(t-1)*(lv2(t-1)-1)*nv2(t-1) - nv2(t-1);

            if pibv2(t) > 0
                rhseev2(t) = exp(x*eta1_rhsee_pol);
                rbv2(t)    = exp(x*eta1_rbar_pol);
                nv2(t) = chi1*( chi0*pibv2(t) + nv2(t-1) ) + n0;
            else
                rhseev2(t) = exp(x*eta3_rhsee_pol);
                rbv2(t)    = exp(x*eta3_rbar_pol);
                nv2(t) = chi1*(      pibv2(t) + nv2(t-1) ) + n0;
            end

        else
            pibv2(t)  = ( exp(rkhat2(t)) + 1-delta )*lv2(t-1)*nv2(t-1) - rbv2(t-1)*(1+lambda*(1-gamma))*(lv2(t-1)-1)*nv2(t-1) - nv2(t-1);
            nv2(t) = n_ss*(1 - nd/100);
            rhseev2(t) = exp(x*eta2_rhsee_pol);
            rbv2(t)    = exp(x*eta2_rbar_pol);
            runv2(t) = 1;
        end

        cv2(t) = rhseev2(t)^(-1) + psi*hv2(t)^(1+1/nu)/(1+1/nu);
        iv2(t) = yv2(t) - cv2(t);
        kv2(t) = (1-delta)*kv2(t-1) + iv2(t);

        lv2(t)  = kv2(t)/nv2(t);

        erk_next2 = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t));

        % check policy
        if lv2(t)>thre_l && xiv2(t)>=1

            lv2(t) = thre_l;
            kv2(t) = lv2(t)*nv2(t);
            iv2(t) = kv2(t) - (1-delta)*kv2(t-1);
            cv2(t) = yv2(t) - iv2(t);
            rhseev2(t) = exp( sum(x*eta4_rhsee_pol ) );
            rbv2(t)    = 1/( cv2(t) - psi*hv2(t)^(1+1/nu)/(1+1/nu) ) / rhseev2(t);
            erk_next2  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                       + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t));

        elseif lv2(t)>thre_l && xiv2(t)<1
            
            lv2(t) = thre_l;
            kv2(t) = lv2(t)*nv2(t);
            iv2(t) = kv2(t) - (1-delta)*kv2(t-1);
            cv2(t) = yv2(t) - iv2(t);
            rhseev2(t) = exp( sum(x*eta2_rhsee_pol ) );
            rbv2(t)    = 1/( cv2(t) - psi*hv2(t)^(1+1/nu)/(1+1/nu) ) / rhseev2(t);
            erk_next2  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                       + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t));
            
        else
            rbv2(t) = fzero('func_Rb',rbv2(t),options,lv2(t),erk_next2,lambda,gamma,delta,sigma_rk);
        end

        rkhats_next2 = log ( rbv2(t) * (1-1/lv2(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp2     = (rkhats_next2 - erk_next2) / sigma_rk;
        cdf2 = normcdf(z_temp2);
        prv2(t) = cdf2;

    % fix tau policy %%%%%%%%%%%%%%%%%%
    
        hv3(t)  = ((1-alpha)/psi*exp(av(t))*kv3(t-1)^(alpha))^(1/(alpha+1/nu));
        yv3(t)  = exp(av(t)) * kv3(t-1)^(alpha)*hv3(t)^(1-alpha);

        rkhats3(t) = log ( rbv3(t-1) * (1-1/lv3(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat3(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv3(t-1));
        xiv3(t)  = exp(rkhat3(t))/exp(rkhats3(t));

        K  = [kv3(t-1) rbv3(t-1) nv3(t-1) xiv3(t)];
        x0 = (log(K)).^expmx;
        x  = prod(x0,2)';
        if xiv3(t) >= 1
            pibv3(t)  = ( exp(rkhat3(t)) + 1-delta )*lv3(t-1)*nv3(t-1) - rbv3(t-1)*(lv3(t-1)-1)*nv3(t-1) - nv3(t-1);

            if pibv3(t) > 0
                rhseev3(t) = exp(x*eta1_rhsee_pol3);
                rbv3(t)    = exp(x*eta1_rbar_pol3);
                nv3(t) = chi1*( chi0*pibv3(t) + nv3(t-1) ) + n0;
            else
                rhseev3(t) = exp(x*eta3_rhsee_pol3);
                rbv3(t)    = exp(x*eta3_rbar_pol3);
                nv3(t) = chi1*(      pibv3(t) + nv3(t-1) ) + n0;
            end

        else
            pibv3(t)  = ( exp(rkhat3(t)) + 1-delta )*lv3(t-1)*nv3(t-1) - rbv3(t-1)*(1+lambda*(1-gamma))*(lv3(t-1)-1)*nv3(t-1) - nv3(t-1);
            nv3(t) = n_ss*(1 - nd/100);
            rhseev3(t) = exp(x*eta2_rhsee_pol3);
            rbv3(t)    = exp(x*eta2_rbar_pol3);
            runv3(t) = 1;
        end

        cv3(t) = rhseev3(t)^(-1) + psi*hv3(t)^(1+1/nu)/(1+1/nu);
        iv3(t) = yv3(t) - cv3(t);
        kv3(t) = (1-delta)*kv3(t-1) + iv3(t);

        lv3(t)  = kv3(t)/nv3(t);

        erk_next3 = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv3(t));

        % check policy
        lv_unreg(t) = lv3(t)/(1-tau);
        kv_unreg(t) = lv_unreg(t)*nv3(t);
        erk_next_unreg(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                          + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv_unreg(t));
        rbv3(t) = fzero('func_Rb',rbv3(t),options,lv_unreg(t),erk_next_unreg(t),lambda,gamma,delta,sigma_rk);

        rkhats_next3 = log ( rbv3(t) * (1-1/lv3(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp3     = (rkhats_next3 - erk_next3) / sigma_rk;
        cdf3 = normcdf(z_temp3);
        prv3(t) = cdf3;

    end
    % End of T period simulation

%% plot

close all

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','output')
plot(-t_before:t_after,yv(2:1+t_before+1+t_after)/avg_y_nopol*100-100,'Linewidth',3)
hold on
plot(-t_before:t_after,yv3(2:1+t_before+1+t_after)/avg_y_nopol*100-100,'k-.','Linewidth',3);
hold on
plot(-t_before:t_after,yv2(2:1+t_before+1+t_after)/avg_y_nopol*100-100,'r:','Linewidth',3);
hold on
plot(-t_before:t_after,zeros(1,t_before+t_after+1),'k:','Linewidth',2);
%ylim([-0.05 0.05]);
%yticks(-0.05:0.01:0.05)
set(gca,'FontSize',30);
xlim([-t_before t_after]);
xticks(-t_before:4:t_after)
xtickangle(0)
grid on
ytickformat('percentage')

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','investment')
plot(-t_before:t_after,kv(2:1+t_before+1+t_after)/avg_k_nopol*100-100,'Linewidth',3)
hold on
plot(-t_before:t_after,kv3(2:1+t_before+1+t_after)/avg_k_nopol*100-100,'k-.','Linewidth',3);
hold on
plot(-t_before:t_after,kv2(2:1+t_before+1+t_after)/avg_k_nopol*100-100,'r:','Linewidth',3);
hold on
plot(-t_before:t_after,zeros(1,t_before+t_after+1),'k:','Linewidth',2);
%ylim([-0.05 0.05]);
%yticks(-0.05:0.01:0.05)
set(gca,'FontSize',30);
xlim([-t_before t_after]);
xticks(-t_before:4:t_after)
xtickangle(0)
grid on
ytickformat('percentage')
legend({'No policy','$$\tau=2.5\%$$','$$\bar{L}=15$$'},'Fontsize',18,'NumColumns',3,'Interpreter','latex')

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','net worth')
plot(-t_before:t_after,nv(2:1+t_before+1+t_after)/avg_n_nopol*100-100,'Linewidth',3)
hold on
plot(-t_before:t_after,nv3(2:1+t_before+1+t_after)/avg_n_nopol*100-100,'k-.','Linewidth',3);
hold on
plot(-t_before:t_after,nv2(2:1+t_before+1+t_after)/avg_n_nopol*100-100,'r:','Linewidth',3);
hold on
plot(-t_before:t_after,zeros(1,t_before+t_after+1),'k:','Linewidth',2);
%hold on
%plot(-t_before:t_after,ones(1,t_before+t_after+1)*n_stss_nopol,'k:','Linewidth',3);
%ylim([-0.05 0.05]);
%yticks(-0.05:0.01:0.05)
set(gca,'FontSize',30);
xlim([-t_before t_after]);
xticks(-t_before:4:t_after)
xtickangle(0)
grid on
ytickformat('percentage')

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','leverage')
plot(-t_before:t_after,lv(2:1+t_before+1+t_after),'Linewidth',3)
hold on
plot(-t_before:t_after,lv3(2:1+t_before+1+t_after),'k-.','Linewidth',3);
hold on
plot(-t_before:t_after,lv2(2:1+t_before+1+t_after),'r:','Linewidth',3);
hold on
plot(-t_before:t_after,ones(1,t_before+t_after+1)*thre_l,'k:','Linewidth',2);
%ylim([-0.05 0.05]);
%yticks(-0.05:0.01:0.05)
set(gca,'FontSize',30);
xlim([-t_before t_after]);
xticks(-t_before:4:t_after)
xtickangle(0)
grid on

return
